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First-principles density functional calculations for graphene and defected graphene are used to 
examine when the quasi-2D electrons near the Fermi energy in graphene could be represented by 
massless fermions obeying a Dirac- Weyl (DW) equation. The DW model is found to be inapplicable 
to defected graphene containing even ~3% vacancies or N substitution. However, the DW model 
holds in the presence of weakly adsorbed molecular layers. The possibility of spin-polarized phases 
(SPP) of DW-massless fermions in pure graphene is considered. The exchange energy is evaluated 
from the analytic pair-distribution functions as well as in fc-space. The kinetic energy enhancement 
of the sipn-polarized phase nearly cancels the exchange enhancement, and the correlation energy 
plays a dominant residual role. The correlation energies are estimated via a model four-component 
2D electron fluid whose Coulomb coupling matches that of graphene. While SPPs appear with 
exchange only, the inclusion of correlations suppresses them in ideal graphene. 
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I. INTRODUCTION. 

Graphene and related materials (e.g, nanotubes, 
fullerenes) have become a mine of novel technologies and 
new horizons in physics [lj]. These include cosmological 
models on honeycomb branes, superconductivity on bi- 
partite lattices 3, nanotubes 0, Q, Hubbard models @, 
spin-phase transitions [3], nanostructuresjH and other as- 
pects of strongly correlated electrons 0. The carbon 
atoms in graphene form a quasi-two-dimensional (Q2D) 
honeycomb lattice and contribute one electron per carbon 
to form an unusual 2D electron system (2DES) with a 
massless Fermions obeying a Dirac- Weyl (DW) equation 
near the Fermi pointsjiol. [Til [l2| . The hexagonal Bril- 
louin zone has two inequivalent points K=(l/3, 1/^/3) 
and K'=(— 1/3, 1/^3), in units of 27r/ao, where ao is the 
lattice constant. The simplest tight-binding model with 
nearest-neighbour hopping t is sufficient to describe the 
valence and conduction bands {tt and tt*) near the K, K' 
points, i.e, at the Fermi energy Ep, where the bandgap 
is zero. The Dirac- Weyl 2D electron system (DW-2DES) 
is nominally "half-filled" , with the tt* band unoccupied, 
(see Fig. [1]) and has spin and valley degeneracies, with a 
Berry phase associated with the valley index 12'. 

The above picture assumes a perfect 2-D sheet of C 
atoms in a honeycomb lattice held in place by the er- 
bonding structure of the C framework. In practice, since 
defects are favoured by the entropy term in the free en- 
ergy, some carbon atoms may be missing, forming va- 
cancies; they may also be substituted with other defect 
atoms. The surface itself may be covered with adsorbed 
gases. Hence the nature of the density of states (DOS) 
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near the Fermi energy in systems with vacancies and sub- 
stituted atoms needs to be considered [IH 03]. Removing 
a carbon atom effectively removes four valance electrons 
from the system, and the resulting vacancy may or may 
not lead to a magnetic, conducting or insulating ground 
state. Replacing a C atom by, say, a nitrogen atom pro- 
vides five valance electrons. We find that at typical con- 
centrations of 3% or more, the Dirac- Weyl picture fails, 
and the Fermi energy moves into a bandgap or to regions 
with a high density of states. However, if less extreme sit- 
uations are considered (e.g, adsorbed gases on graphene), 
the DW picture is found to hold true. 

The Fermi liquid found in metals and semi-conductor 
interfaces is characterized by the Wigner-Seitz radius r s 
of the sphere (or disk in 2D) which contains one electron. 
When expressed in atomic units (involving the effective 
mass m* of the electrons), r s becomes the ratio of the 
potential energy to the kinetic energy, i.e, a measure of 
the strength of the Coulomb interaction. Hence r s < 1 
provides a regime where Fermi-liquid perturbation theory 
is strictly valid. The vanishing of the density of states and 
the effective mass of the electrons in DW-2DES near the 
Fermi points implies that the strict Fermi-liquid picture 
breaks down in perfect graphene. 

The DW-2DES can be doped to contain carrier elec- 
trons (in the tt* band), or holes as well (in the tt band), 
and provides a rich system which retains the massless- 
fermion character, as long as the material is not strongly 
modified. However, electron-electron interactions in the 
DW-2DES may modify the tt and tt* bands, lift the sub- 
lattice (valley) degeneracies, or stabilize spin-polarized 
phases (SPP) in preference to the unpolarized state, if 
the doping level or the strength of the Coulomb inter- 
action could be varied. The effect of electron-electron 
interactions has been most extensively studied in Fermi- 
liquid- like electron systems found in GaAs/AlAs inter- 
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FIG. 1: (Color online)Linear dispersion bands near a K point 
where the n* and 7r bands cross. In (a) we show a doped 
unpolarized system with equal occupation of the up-spin (a u ) 
and down-spin states. In (b) the polarized system has only 
electron carriers. In (c) both electron and hole carries occur. 
This is the only possibility for spin polarization if doping is 
zero (&f = and k u = kd)- 



faces or Si/Si02-inversion layers. The SPP in such 
GaAs/Al As-based 2DES, predicted to appear at low cou- 
pling (r s ~ 2 — 4) if perturbation methods are used, get 
pushed to high coupling (?% ~ 26) if non-perturbative 
theory (NPT) were used [H Ezj]. The two- valley 2DES 
does not show a SPP in NPT, unlike for one-valley sys- 
tems, presumably because of the preponderance of the 
(three times as many) direct Coulomb interactions over 
the exchange interactions [Hj]. The exchange and cor- 
relation energy E xc in the 4-component Si/Si02 elec- 
tron system were calculated in Reflla using the classical- 
map hyper-netted-chain (CHNC) technique, accurately 
recovering the Quantum Monte-Carlo (QMC) results in 
every coupling regimes 19]. CHNC provides the pair- 
distribution functions (PDFs) function of the 
coupling strength. Then E xc is evaluated via a cou- 
pling constant integration, providing a fully non-local, 
transparent approach. The method has been success- 
fully applied to the 2DES|p, 3DES[2l[, the 2- valley 
2DES in SiC-2 interfaces [1^, the thick quasi-2DES in 
HIGFETS[13], and the electron- proton system [H]. How- 
ever, the full non-local treatment of exchange and cor- 
relation in graphene involves an 8 x 8 matrix of two- 
component PDFs because of the spin, and valley indices 
as well as the presence of 71", 7r* bands. Hence in this 
study we first consider the exchange energies via an an- 



alytic evaluation of the non-interacting PDFs. The cor- 
relation energies are estimated by appealing to our re- 
sults for the spin-polarized four-component two- valley 2- 
D electron system (2v-2DES) of Si-MOSFETS. The non- 
interacting PDFs of the DW-2DES, G°j(r) involve two 
components, the first being a Bessel function as in the or- 
dinary 2DES, while a second, associated with the cosine 
of the angle of e-e scattering, involves products of Bessel 
and Struve functions. The exchange energy enhancement 
of the sipn-polarized phase nearly cancels the kinetic en- 
ergy enhancement, implying that the correlation energy 
plays the dominant residual role. We find that while 
there are stable SPPs in an exchange-only approach, in- 
cluding the correlation energy using the 2v-2DES data 
stabilizes the graphenc-2DES in the unpolarized state. 
This conclusion is not surprising since the Coulomb cou- 
pling strength in graphene is ~ 2-3, and no SSPS are 
found in electron liquids at such low coupling, except as 
artifacts of low-order theories. 

In the following section we present first-principles cal- 
culations for pure graphene, graphene with ~ 3%, and 
~ 12 % vacancy concentrations, and show that the DW- 
model is inapplicable to such systems. We also consider 
systems with N-substitution as the lattice distortion ef- 
fects are smaller here. Nevertheless, even here the DW- 
model does not seem to be applicable. However, if we 
consider a graphene layer having a metastable sheath of 
N 2 molecules adsorbed on it, with no disruption of the a- 
bonding network, the DW-model does remain applicable. 
Thus, having established the limits of the DW-model, we 
proceed to examine the exchange and correlation effects 
among DW-fermions, and show that a stable spin phase 
transition is found only in "exchange only" models which 
neglect electron correlation effects. 



II. DENSITY-FUNCTIONAL CALCULATIONS 
OF GRAPHENE SYSTEMS 

Simple tight-binding models (TBM) can be trusted 
only if they are validated by more detailed calcula- 
tions and experiments. While TBM can be success- 
fully exploited within a limited energy window for pure 
graphene, the effect of vacancies and lattice substituants 
etc., needs detailed consideration. A vacancy removes 4 
valence electrons, distorts the bond lengths and angles 
around the vacancy, creating pentagonal networks and 
compensating larger networks, localizing electrons near 
the vacancy and changing the structure of the electronic 
density of states (DOS). The bond- lengths and the net- 
work structure is better preserved with N substitution. 
Here an extra electron is added to the graphene system 
for every N substitution. Attempting to treat such ef- 
fects using tight-binding methods augmented by, say, T- 
matrix theory etc., to account for impurity effects are 
well known to be strongly model dependent. Neverthe- 
less, some insight has already been gained from short- 
ranged scatterer models where lattice relaxation and 
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other very important issues are not handled [3l|, [32], |33j 
However, density-functional theory (DFT) has an excel- 
lent track record in just such problems where electronic 
and ionic energy minimization can be carried out until 
the Hellman-Feynman forces on atoms around the va- 
cancy or the substituant are reduced to zero. There 
are a number of such DFT calculations already in the 
literature [ID, [3], mainly concerned with energetics and 
bonding. Here we examine the bands and DOS of these 
systems, with an eye on the limits of validity of the DW- 
2DES model. 

We have used the Vienna ab initio simulation pack- 
age (VASP)[25| which implements a spin-density func- 
tional periodic plane-wave basis calculation. The pro- 
jected augmented wave (PAW) pseudopotentials[25[ have 
been used for C and N. The C pseudopotential has al- 
ready been used in several graphene-type calculations 
(e.g., Ref. fl3h . The N pseudopotential was also further 
tested by a study of the N2 dimcr where an equilibrium 
bond distance of 1.11 Angstroms was obtained, in good 
agreement with other DFT calculations as well as results 
from detailed configuration interaction studies etc.[26|. 
where a value of 1.095 A is reached. The ~ 3% and ~ 
12% vacancy calculations were done with 32-atom and 8- 
atom graphene-like unit cells. These systems are thus not 
truly disordered, but provide reference densities of states 
which acquire smearing when some disorder is introduced 
by s amp ling other configurations using larger simulation 
cells[23|. When vacancies are introduced into graphene, 
the large stresses are relieved by the neighbours (C-atoms 
near the vacancy) moving towards the vacancy. The 
distortion persists to at least the third neighbours, and 
generates a bond-length distribution varying from about 
1.398 A to 1.45 A. If the vacancies are replaced by N 
substitution, the structural distortions are smaller but 
the changes in the DOS can be equally drastic, as we 
show below. 

Yuchen Ma et al.[l3| have found that N2 molecules 
may form a metastable layer near carbon nanotubes and 
graphene. Iyakutti et al[27J have also found similar sta- 
bilization, where the N2 layer is less stable than if it were 
at infinity, but held in place by an energy barrier. Al- 
though this is a 'weakly adsorbed' state, the interaction 
energy between the graphene layer and the N2 layer is 
about as strong as between two graphene sheets. Hence 
we have looked at the band-structure and DOS of such 
a fully N 2 covered graphene sheet as well. Here the 
Dirac-Weyl behaviour is preserved. In Fig. [5] we show 
the band structure of pure graphene (top panel), and 
also a graphene sheet with a layer of N2 molecules, with 
each N2 aligned on every Kekule-bond position (bottom 
panel). The sheet of N2 molecules is positioned about 
3 A above the graphene sheet. The interaction energy 
per carbon (or per N) is about 0.2 eV. The regime of D- 
W linear dispersion around the K point is reduced from 
that of pure graphene. This is in contrast to the interac- 
tion between two graphene sheets, where the bands near 
the K-point become parabolic [l5j. 
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FIG. 2: (Color online) Band structure of pure graphene (top 
panel), and that of fully N2 covered graphene. The Fermi 
energy is at E+i 7 =0. Conduction bands are shown in red. 
The bands obtained from a standard sp 3 tight-binding scheme 
are compared with the DFT bands for graphene in the top 
panel. 



However, when vacancies or N atoms are introduced 
into the graphene network at relatively significant con- 
centrations (e.g., 3%-12 %), the DOS modifies and shifts 
to produce a large density of states near the Fermi en- 
ergy Ep. This is shown in the top panel of Fig. [3] The 
calculation is for an ordered system with 1/8 or 1/32 al- 
loying of vacancies or N, and may be used as reference 
points for a consideration of a disordered system. The 
calculations used exchange-correlation functionals based 
on the Ceper ley- Alder type: 29], and no spin polarized 
states are expected for a conductive system with an en- 
hance DOS, as seen in the figure. On the other hand, the 
existence of such an enhanced DOS near Ep in graphene 
with 12% vacancies, would be of importance to possible 
superconductivity of these systems. 

The electronic density of states for a lower density 
vacancy system, i.e., with one vacancy per 32 carbon 
atoms, is shown in the lower panel. Here the system ac- 
quires a gap near Ep, and the material is an insulator. 
If the vacancy concentration is further decreased, we ex- 
pect the energy gap to slowly decrease and recover the 
Dirac-Weyl model in the limit of pure graphene. But, 
as seen from the DOS at 12% vacancies, we see that as 
the vacancy concentration is increased from 3%, the en- 
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ergy gap closes and Ep positions to a high-DOS region. 
Meanwhile, an energy gap opens about 1 eV below the 
Fermi energy. Thus we see that a vacancy induced metal- 
insulator transition is possible (see however, Ref. l30h - 
This picture becomes considerably less sharp if the va- 
cancies are not considered to form a periodic array. In 
any case, our conclusion that the Dirac-Weyl model is 
inapplicable even at 3% vacancy concentrations proba- 
bly remains valid. The observations of the quantum Hall 
effect and other signatures of the DW model clearly in- 
dicate that good low-vacancy regions of graphene foils 
are the subject of these experiments. Another aspect of 
low-density vacancies or substitucnts in graphene is the 
issue of spin-polarized ground states. The reported re- 
sults depend on the size and edge structure of finite sheet 
fragments (HI, or the possibility of further C-adatom ad- 
sorption on N-substituted sites [HI]. The calculations are 
sensitive to the treatment of exchange and correlation, 
as in first-principles theories of magnetism in transition- 
metal oxides. Hence a reliable discussion of magnetism 
in graphene-like systems would require an assessment of 
Coulomb interactions in Dirac-Weyl electrons. Hence 
we revert to the main object of the present study, viz., 
the nature of exchange and correlation effects in pure 
graphene. 
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III. THE 2DES IN GRAPHENE NEAR THE 
FERMI POINTS. 

The band structure of graphene near the K points, i.e., 
close to the Fermi energy, are displayed in Fig. [2] The 
kinetic energy near the K points is given by a Dirac-Weyl 
Hamiltonian of the form: 



H k = V F (p x T z a x +p v o-y) 



(1) 



Here r 2 = ±1 defines the degenerate valleys, and o~ x , a y 
denote the x and y Pauli matrices that act in the space of 
the two atoms in each unit cell. The 7T, 7r* bands of spin 
and valley degenerate states (Fig. [T]) show a linear disper- 
sion E — ±Vphk. This form requires a cutoff momentum 
K c such that the number of states in the Brillouin zone 
is conserved. That is, if Aq is the area per carbon, then 



K 2 C = 47r(l/A ) 



(2) 



. The electron density iV c at half-filling is 1/Ao, with 
Aq = Oq^/3/2, since one ir electron of arbitrary spin is 
provided by each carbon atom. The Fermi velocity Vp = 
tao^3/2 is thus the slope of the linear dispersion, with 
V F ~ 5.5 eVA. If the DW-2DES is embedded in a medium 
with dielectric constant eg, then we define 
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e 2 /e 

hV F 



eoao 



7(V3/2). 



(3) 



This is the ratio of a typical Coulomb energy to the hop- 
ping energy and hence is usually taken as the Coulomb 
coupling constant of the DW-2DES. This ratio plays the 



FIG. 3: (Color online)DOS of pure graphene compared with 
~ 3% and ~12% concentrations of substituted N atoms(top 
panel), calculated using a periodic substitution model. The 
Fermi energy is set to zero. Bottom panel shows the effect of 
~ 3%, 12% vacancies in graphene. 



same role as the r s parameter in electron-gas theory of 
nonrclativistic finite-mass fermions, and is a measure of 
the strength of the Coulomb interaction. The usual r s is 
not available for DW-2DES since the effective mass m* is 
zero and there is no effective Bohr radius. The coupling 
constant g° is maximized if eo is unity, and consistent 
with this case we assume g° = 2.7, Vp = 5.39 eVA , with 
e 2 /e =14.4, for our DW-2DES studies. 

The 4-component-eigenfunction envelopes of the ki- 
netic energy term are made up of two-component func- 
tions U = (6,e^ fe ), U' = (e^<" b ) and O = (0,0) where 
4>k is the angle of the vector k in the 2-D plane. Thus 



Ff ? (r) = (2A)^(U,0) TXa 
{2Afl\0,U') TX o 



b,k 



(4) 
(5) 



Here b = ±1 is a it* , ir band index, (• • • )t indicates the 
transpose, and \a is the spin function. Then, using v = 
1,2 as a valley index, the Coulomb interaction term in 



5 



the Hamiltonian may be written in the form: 



1 

8A 



E E V " [&i6 4 e* w * (k) -^ k+q)} + l 

Vi ,bi ,<Ti k.p.q 

x ^ 2 & 3 e i{0 * (p) ~ 0(p+q)} + 1 x 

(6) 

Here a + , a are electron creation and annihilation opera- 
tors and V q = 2ire 2 / (e^q) is the 2D Coulomb interaction. 
The phase factors introduce a novel cos (9) contribution 
where 9 is the scattering angle, not found in the usual 
jcllium-2DES. The resulting form of the exchange energy 
per Carbon is: 



E x /E u 



A g°/k c 1 
(2tt) 2 4 



xkp 



E 

l + 6i6 2 cos(0) 



2ti 



|k-p| 



d9dkdp 



n burT (k)n b2iCr (p) 



(7) 



In the above we have included the intrinsic coupling con- 
stant g° and the energy unit E u — Vpk c in the expres- 
sions. Here k c — K c / \/2 = y/(4irn c ) is based on the 
electron density per spin species, n c — N c /2 = 1/(2 An). 
The above form of the exchange energy can be reduced 
to an evaluation of a few elliptic integrals |7j. The nor- 
mal "half-filled" DW-2DES can be doped with electrons 
or holes; but it is easy to show that symmetry enables 
us to limit to one type of doping. However, given a sys- 
tem, with an areal density of N$ dopant electrons per 
valley, with ng — N$/2 per spin, the carriers in the spin- 
polarized system could be electrons only, or both elec- 
trons and holes, as shown in Fig. Q] for the tv* and tv 
bands at one K point. The intrinsic system with ns = 
can be an unpolarized state, or spin-polarized state with 
electrons and holes. Such purely exchange driven systems 
have been studied by Peres et al. 0] , while the correlations 
effects have not been considered. Here we evaluate the 
exchange energy E x from from the non-interacting PDFs, 
and include the correlation energy E c estimated from the 



2v-2DES with the same coupling strength (r s 
spin-polarization. 



and 



A. Electron Pair-distribution functions of the 
Dirac-Weyl 2DES 

Although we are dealing with an intrinsically four- 
component system (2- valleys, 2-spin states), as seen from 
Fig. [TJ we need to consider the redistribution of electrons 
and holes among the tv* and tv bands when comparing the 
energy of spin-polarized states with the corresponding 
unpolarized state. The exchange energy is a consequence 
of the anti-symmetry of the Slater determinant made up 
of non-interacting eigenfunctions. Thus only the non- 
interacting PDFs are needed to calculate the exchange 
energy. They are also the spring-board for calculating 
the interacting PDFS via the CHNC method. 



BO 
II 





1 


\\ / / 

\ \ / / / 

•% \ / / 




/ intraband 


h u (r) 


/ intraband 


h>) 


/ interband 


h S 13 (r) . 


/ spin-parallel PDFs. 
/ deviation from 1/2 filling, 


- 

n s /n =0.2 

be _ 


/ zero spin polarization. 

y 


1 



r/a n 



FIG. 4: (Color online)The Bessel-like and Struve-like non- 
interacting, parallel-spin PCFs h b (r) and h"(r) for the unpo- 
larized doped system. The bands are numbers as in Fig. [Hb) . 
The anti-parallel non-interacting PCFs are zero. The lattice 
constant ao = 2.47A. 



From Fig. Q]we see that the e-e interactions at a given 
valley can be constructed from (i) interactions with a 
tt*(t u band of up-spin electrons of density n u , filled to 
k u , (ii) a n* ad set of electrons or a ivaa spin-down holes, 
of density n<j, filled up to kd (iii)the iva u band, with 
electron density n c , filled to k c , and (iv)the nad band, 
density n c , filled to k c There will also be similar inter- 
valley terms. Each term in this 4x4 matrix, denoted 
by Qij(r) where i, j = 1 ■ ■ -4, will have two components 
associated with those in U' and U = (b,e l ^ k ). Thus 
Gij( r ) — 9ij ( r ) > 9tj ( r ) j where the superfixes 'bjs' indi- 
cate that the noninteracting forms are Bessel-function 
like, and Struve-function like, respectively. These com- 
ponents will be denoted by a superfix c = b, s. The 
Struve form arises from the cos(0) terms in the Coulomb 
interaction. The numbering scheme of the matrix is 
shown in Fig. [ljb). We denote the pair-correlation func- 
tions (PCFs) Hij(r) = Qij(r) — 1, or the components 
by Mj( r ) ~ 9ij( r ) ~ 1- The non-interacting forms are 
indicated by a superscript zero. Thus we have: 



j 0,6/ \ 

K ( r ) 



= -(UiUj) 

2 



dkx f* dk 2 <(ki 



(ki-k 2 )-: 



o (2TvfJ (2tvY 
2 



- — Ji(k t r)- — Ji(kjf) 
k,r k,r 



(8) 



h °if( r ) = -(niUj)- 



dk x [*» dk 2 



cos(fi — f 2 Je 

TV TV 



o (2TvfJ (2tvY 

(k!-k 2 )-r 



fcjr kjr 



[JqH\ — JiH ]i[JoHi — H n Ji]j 
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Here Jo, Ji are Bessel functions, while Hq and Hi are 
Struve -functions. Also, in [J Hi — JiH } i the functions 
are evaluated at the argument for. That is, 

[Jo-Hi - JiHoji = JoirkjHxirki) - Ji(rfo)ff (rfo) (9) 

The wavevectors fo = y / (47rnj) for each component i, of 
density Uj. We show typical noninteracting PCFs for a 
doped, unpolarized case as in, Fig.JTJa), with the doping 
fraction ng/n c = 0.2. In CHNC, the exchange-hole is 
mapped exactly into a classical Coulomb fluid using the 
Lado procedure [21]. The figure shows that the exchange- 
hole is strongly reduced by the presence of the cos(#) 
term which has been averaged into the Struve-like PCFs 
h s (r). When the Coulomb interaction is included, the 
cos(#) term has a similar mitigating effect and exchange- 
correlation in the DW-2DES is considerably weaker than 
in the corresponding 2-valley 2DES. The CHNC calcu- 
lation for the 2v-2DES for the conditions stipulated in 
Fig. 2] show that the correlation energy is only about a 
third of the exchange energy. This motivates our use of 
the 2v-2DES for the correlation energy, while the E x is 
exactly evaluated. 

B. The kinetic and exchange energies. 

When the doping per valley is N$ — 2ns, the total 
number of electrons per valley is N t = N c + N$. Also, us- 
ing the i — 1, 2, 3, 4 notation of Fig.QJb), we set n\ = n u , 
n 2 = n&, U3 = 71.4 = n c . Hence the spin polariza- 
tion s — n u — b d n d , where the band index bd = — 1 
for holes. The degree of spin-polarization ( — s/N t . 
The composition fractions, inclusive of the valley index 
v=l,2 are x V i = n,i/2N t . We note that k F — y/(2wns), 
k u = y/{2n(ns + s)}, kd = y/{2ir\ns — s\}. The exchange 
energy E x (ng, C) can be written as: 

— ^— ( 10 ) 

ij 

It is understood that the Struve-like component in Q v ,v,ij 
where v labels the valleys, is summed with the appropri- 
ate bibj band ± factors. Only a sum over the components 
in one valley is needed in evaluating the exchange. The 
above formula can be made more explicit, by introduc- 
ing the b^j band ± factors in each case. Thus the total 
kinetic and exchange energy Eu x = KJE+E X , for the case 
(b) of Fig. [T] can be written in terms of n F ,n u ,nd and n c 
as in Eq. 1101 or in terms of kp,k u ,kd, k c and Aq as: 

E kx (Q = ^hV F (k 3 u + k 3 d )- (li) 

D7T 

^( 5 ^y F /4)(7r/2)[fc^H(r) + k\U 22 {r) + 
2k 2 c {k 2 u n 13 (r) + k 2 d H 14 (r)}} 



The kinetic and exchange energy, E tx (C = 0), of the 
unpolarized system, shown in Fig. [TJa), is obtained by 
setting fcf = in Eq. [TTJ The exchange energy for the 
case involving both electron and hole carriers, Fig.[TJc) is 
given by a small change in Eq. 111! where the last term 
changes sign and becomes — k^Hi^r) . If this is simpli- 
fied and written using elliptic integrals, as in Ref. 0, the 
change in the exchange energy for electron and hole car- 
riers, compared to the unpolarized case becomes 

E X (C) = E x {n 5 ,C = 0)+AE x (C) (12) 

A£UC) = j^(g /^(Eu/k c )[(k 3 u + k 3 d -2k F )R 1 (l) + 

2k c k 2 u R 2 (k u /k c ) ~ 2k c k 2 d R 2 {k d /k c ) 
-Ak c k 2 F R 2 {k F /kc)] (13) 

As before, E u = WJpk c is the unit of energy. For com- 
pleteness, we note that 

E x (n s = 0,C = 0) = - 7 ^(g a /2)E u k 2 c Ri(l) (14) 
(2tt) 2 

R X {1) = 3.776 (15) 

Here the energy is per carbon atom and we have used the 
notation R\ (x) , R 2 (x) for the elliptic integrals as given 
m Ref. 0, and taken the opportunity to correct a few 
typograpical errors in their Eq. (15). The energy differ- 
ence which determines the competition between the un- 
polarized and polarized phases, i.e., AEk x (C)i includes 
the kinetic-energy corrections as well as the change in 
the exchange energy. It is plotted in Fig. [5] We have 
done the calculations in r-space using the PDFS, and in 
fc-space via elliptic integrals, to provide independent nu- 
merical checks. Fig. [5] shows that stable spin-polarized 
phases appear in electron-carrier systems, if kinetic and 
exchange energies are used in the total energy. A note- 
worthy feature of Fig. [5] is the strong cancellation of the 
kinetic energy by the exchange energy, leading to net en- 
ergies which are about 1% of the energy scale E u = Vpk c . 
Thus the stage is set for the phase stabilities to be de- 
termined by the correlation energies which are left out in 
Fig. E 

IV. CALCULATIONS INCLUDING THE 
CORRELATION ENERGY. 

A rigorous calculation of the correlation energy re- 
quires the self-consistent evaluation of 8 x 8 matrix of 
pair-distribution functions for many values of the cou- 
pling constant < A < 1 in the interaction A/r, and an 
integration over the coupling constant A. 

Although this can be envisaged within the CHNC ap- 
proach, it still remains a very arduous task. Even when 
all the inherent symmetries in the problem are taken into 
account, some two-dozen PDFs need to be evaluated. In- 
stead we outline a simplified scheme where we use the 
results for the correlation energy of the two-valley, two- 
spin (i.e., 4 component) jellium 2DES to reconstruct the 
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FIG. 5: (Color online)Upper panel- the energy difference 
AEkx, i.e., K.E+exchange, between the polarized and un- 
polarized phases, in units of E u = Vpfc c , as a function of 
the spin polarization £ and the dopent density rid- Electron- 
carrier systems, Fig. QJb) are more stable than electron-hole 
systems, and show stable spin-polarized states. However, ad- 
dition of the correlation energy (see Fig. [6j makes the unpo- 
larized state the most stable phase. 



Spin polarization (Q 
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FIG. 6: (Color online)Upper panel- The correlation-energy 
difference AE c (rid, £), i.e., between the polarized and unpo- 
larized phases, in units of E u = Vf/Cc, as a function of the 
spin polarization £ and the dopent density rid- The lower 
panel shows the total energy difference between the polarized 
and unpolarized phases. The unpolarized state is the most 
stable phase. 



correlation energy of the DW-2DES for equivalent values 
of the ratio of the Coulomb interaction and the kinetic 
energy, as discussed below. The method we use here 
provides a different, possibly more transparent approach 
than that given in a previous discussion [34|. The conclu- 
sions from the two methods are in agreement. 

When the electron density of the jellium 2v-2DES is 
changed, the r s value changes. In contrast, when the elec- 
tron density (or, equivalently, the number of electrons per 
carbon) in the DW-2DES of perfect graphene is changed, 
the coupling constant g° remains unchanged. Hence the 
correlation energy per carbon for a doping situation in- 
volving a total of tit electrons per carbon, at a spin- 
polaraization ( might be approximated by riTe 2 c v {r Sl Q, 
with r s — g° , where e 2v (r s ^) is the corresponding 2v- 
2DES correlation energy. However, this result is in terms 
of an effective atomic unit E m , intrinsic to the jellium- 
2DES, such that the Fermi energy Ep — E eu /r 2 . Since 
we identify the coupling constant g° to be r s , the effec- 
tive atomic unit E au = Epr 2 . A full evaluation of the 
4-component jellium 2DES from the interacting PDFs 
has been carried out and the correlation energy is given 
in parametrized form in Eq. 5 of Ref. [l8l . In transfer- 
ring from the 2v-2DES to the DW-2DES we note that 
e 2 /eo which is unity in 2v-2DES becomes <?°Vf in the 
DW-2DES. The correlation energy enhancement is: 

AE c (n d ,0 = E c (n d ,()- E c {n d ,0) (16) 

The correlation energies in jellium-2DES are evaluated 
from PDFs whose non-interacting forms are Bessel-type 
PDFs, where as the DW-2DES contain only about 1/2 
the number of Bessel-type PDFS, while the other are 
Struve-type PDFs which bringing in a minor contribu- 
tion. Thus an upper bound would be to use the estimate 
of AE c (n d , C) given above, while a lower-bound would be 
about 1/2 the above estimate. The AE C calculated from 
the jellium 4-component system using the above scheme 
is shown in Fig. [5] The total energy difference inclusive of 
the kinetic energy, E x , and the estimated E c , between the 
polarized and unpolarized phases is shown in in the lower 
panel of Fig. Sice the kinetic-energy enhancement of 
the polarized phase is compensated by the exchange en- 
hancement, the total energy enhancement is determined 
by the correlation effects. Thus we see that the inclusion 
of the correlation energy suppresses the spin-polarized 
phase found in the exchange-only calculation. 

In this work we have kept the Coulomb coupling fixed 
at g° = 2.7 typical of graphene, unlike in other studies^, 
0] where the coupling strength g is taken as a tunable 
parameter, in the spirit of Hubbard-model studies. If 
dielectric screening is taken into account (3o| . the coupling 
is reduced and many-body effects become smaller. We do 
not see a practical experimental scheme for increasing the 
value of the Coulomb coupling strength g° beyond 2.7 in 
the graphene system. 

Even in the one- valley 2DES, the SPP of low-order the- 
ories is pushed to g ~ 26 — 27. In the 2v-2DES, direct 
terms predominate over exchange interactions, and the 
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SPP is not found in CHNCpji] or calculations. graphene-2DES as well. This result is in agreement with 

We see that the inclusion of correlations within a reason- the conclusions of Ref. [li- 
able scheme suppresses the exchange-driven SPP in the 
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